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Abstract 

In this paper, we study a strongly coupled two-prey one-predator system. We first prove 
the unique positive equilibrium solution is globally asymptotically stable for the cor¬ 
responding kinetic system (the system without diffusion) and remains locally linearly 
stable for the reaction-diffusion system without cross-diffusion, hence it does not belong 
to the classical Turing instability scheme. Moreover we prove that the positive equilib¬ 
rium solution is globally asymptotically stable for the reaction-diffusion system without 
cross-diffusion. But it becomes linear unstable only when cross-diffusion also plays a role 
in the react ion-diffusion system, thus it is a cross-diffusion induced instability. Finally, 
the corresponding numerical simulations are also demonstrated and we obtain the spatial 
patterns. 
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1. Introduction 


In 2009, Elettreby considered the following prey-predator model |3j 
r 

u\ = aui(l - Mi) - M1M3 := Mi/i(wi, m 3 ), 

< u ’ 2 = bu 2 (l - u 2 ) - u 2 u 3 := m 2 / 2 (m 2 ,m 3 ), (1-1) 

M3 = -CM3 + ( du i + eu 2 )u 3 ■■= u 3 f 3 (ui, m 2 , m 3 ), 
where Mi,m 2 and u 3 are the population densities of three species. This system models the 
dynamic of two-prey one-predator ecosystem, i.e. the third species preys on the second 
and the first one. In the absence of any predation, each term of preys grows logistically. 
The effect of the predation is to reduce the prey growth rate. In the absence of any prey 
for sustenance, the predator’s death rate results in inverse decay, which is the term — cm 3 . 
The prey’s contribution to growth rate of the predators are respectively du±u 3 and cm 2 m 3 . 
They studyed the global stability and persistence of the model. 

However, in reality, individual organisms are distributed in space. We can use the 
reaction-diffusion equations to establish spatio-temporal dynamical system which can 
model the pursuit-evasion phenomenon (predators pursuing prey and prey escaping preda¬ 
tors) in the prey-predator system. Therefore, in present paper we further investigate the 
following reaction-diffusion model with cross-diffusion: 


M U - Apn + fci 3 M 3 )Mi] = aMi(l - Mi) - M 1 M 3 , 

u 2t - A [{k 22 + k 23 u 3 )u 2 ] = bu 2 { 1 - m 2 ) - m 2 m 3 , 

< u 3t - A[(/c 3 iMi + k 32 u 2 + k 33 )u 3 ] = cm| + ( du i + eu 2 )u 3 , 

dui du2 dll3 n 

dr/ dr/ dr/ ’ 

Mi(x,0) = U lo (x),U 2 (x,0) = U 20 (x), U 3 (x, 0) = M3o(x), 


in f 2 x ( 0 , 00 ), 
in Q x ( 0 , 00 ), 
in x ( 0 , 00 ), 
on dfl x ( 0 , 00 ), 

in 


( 1 . 2 ) 


where is a bounded domain in R N with smooth boundary dfl. r) is the unit outward 
normal to dtt. The homogeneous Neumann boundary condition indicates that there is 
zero population flux across the boundary. The parameters a, 6 , c, d, e and k i:J (1 < i, 
j < 3) are all positive constants. k tl is the diffusion rate of f-th species. This diffusion 
term represent simple Brownian type motion of particle dispersal. k tJ (i 7 ^ j) is the cross¬ 
diffusion rate of i-th species. It is necessary to note that the cross-diffusion coefficient 
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may be positive or negative. The positive cross-diffusion coefficient represents that one 
species tends to move in the direction of lower concentration of another species. On the 
contrary, the negative cross-diffusion coefficient denotes the population flux of one species 
in the direction of higher concentration of another species. Here the cross-diffusion term 
presents the tendency of predators to avoid the group defense by a large number of prey, 
i.e. the predator diffuses in the direction of lower concentration of the prey species. More 
biological background can be found in l|, [l4|, |l7|. 


As we know, the problem of cross-diffusion was proposed first 


jy Kerner J 9 ] and first 

applied to competitive population systems by Shigesada et al. [20j. Since then, the role 


of cross-diffusion in the models of many physical, chemical and biological processes has 
been extensively studied. In the field of population dynamics some models of multispecies 
population are described by reaction-diffusion systems. Jorne {§] examined the effect of 
cross diffusion on the diffusive Lotka-Volterra system. They found that the cross-diffusion 
may give rise to instability in the system, although this situation seems quite rare from an 
ecological point of view. Gurtin |(3] developed some mathematical models for population 
dynamics with the inclusion of cross-diffusion as well as self-diffusion and showed that the 
effect of cross-diffusion may give rise to the segregation of two species. Some conditions for 


the existence of global solutions have been given 
[2], Kim 


)y several authors, for example, Deuring 


IP], Pozio and Tesei 19], Yarnada 25]. Moreover, due to a most interesting 


qualitative feature: pattern formation induced by cross-diffusion effect there are some 


works on the diffusion driven instability (Turing instabili 
a non-constant stationary solution, please refer to jll, 12 
references cited therein. 


Y 


22 [) and 


13 


, m Q, 


he existence of 


18 


23 ] and the 


The main purpose of this paper is to study the Turing instability which is driven 
solely from the effect of cross-diffusion by using mathematical analysis and numerical 
simulations. The rest of this paper is organized as follows. In section 2 we show that the 
unique positive equilibrium of the ODE system (11.11) is globally asymptotically stable. In 
the section 3 we show that the positive equilibrium remains linearly stable in the presence 
of self-diffusion. It becomes linearly unstable with the inclusion of some appropriate cross- 
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diffusion influences. The Turing instability occurs only when the cross-diffusion rates k 23 
and k 32 are large. The resulting patterns are computed by a numerical method and also 
we devoted to some conclusions in section 4. 


2. Stability of the positive equilibrium solution of the ODE system 

In this section, we consider the stability of the positive equilibrium solution of the 
system (inj. It is easy to know that if 


abc > max{e(6 — a), d(a — b)} 


(2.3) 


the ODE system (II. ip has a unique positive equilibrium u = (mi, m 2 , M3) which is given by 


u 1 


abc + ae — be _ abc + bd — ad _ ab(d + e) 

-, u 2 =-, M3 =---. 

abc + bd + ae abc + bd + ae abc + bd+ ae 


(2.4) 


We have the following result: 


Theorem 2.1. The unique positive equilibrium u is globally asymptotically stable for the 
ODE system (11.11) . 

Proof. In order to prove the theorem, we need construct a Lyapunov function for the 
system ( 11 . 11 ) . 

V (u (t)) = d(u\ - Ml - Ml In —) + e(n 2 - u 2 - u 2 In —) + (m 3 - m 3 - m 3 In —). ( 2 . 5 ) 

Mi M 2 M 3 

= 0 and E(u) > 0 if u 7 ^ u. By using (II.ip . we compute 
ci(l- )u l + e(l- )m 2 + (1-)m 3 

Ml M 2 M 3 

d(u\ - Mi)[-a(Mi - Ml) - (m 3 - m 3 )] + e(u 2 - u 2 )[-b(u 2 - u 2 ) - (m 3 - m 3 )] 

+ (m 3 - m 3 )[-c(m 3 - m 3 ) + d(ui - Mi) + e(M 2 - m 2 )] 

—ad(ui — Mi) 2 — be(u 2 — u 2 ) 2 — c(m3 — M3) 2 < 0 

for all u 7 ^ u. By the Lyapunov-LaSalle invariance principle Q, u given by (j2,4[) is 
globally asymptotically stable for the kinetic system ( 11 . 11 ) . 


Then E(u) 

dV_ 
d t 
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Theorem 2.2. The unique positive equilibrium u is globally asymptotically stable for the 
reaction-diffusion system (11.21) without cross-diffusion, i.e. kij = 0 for i j. 


Proof. To study the global behavior of system (1.2), we introduce the following Lyapunov 
functional 


W(t) — / V(u(x,t))dx 


( 2 . 6 ) 


where E(u(a;, t)) is given by (12.51) . By direct computation, we have 

[ A V dVi A 

nr = y n grad ^'* dl 

= j (d{ 1 - — ),e(l - —), (1 - —)) ■ (kjiAuj + u,f l ,k 2 2Au2 + u 2 f 2 , 

Jn V u i u 2 u 3 J 

k 33 Au 3 + u 3 f 3 )dx 

= J d ^fcn(l - )A«i^ d x + J e ^fc 22 (l - -^)Aw 2 ^ dx 

+ L{ k 33 il ~^ )Au 3 ) Ax+ L'^ dx - 

From Green’s identity, it follows that 

ka(l - — )Au i ') dx = [ ka(l - — )^ 1 dS - [ k zi X 7 x (l - —) • V x Uidx 

in \ u i J Jon u i on J ^ u* 

= ~ k ii u i uf 2 \'V x u i \ 2 dx < 0 . 

Jn 

Since yjj < 0, f Q < 0. Thus, ^ < 0 for all u 7 ^ u. By the Lyapunov-LaSalle 
invariance principle [7], u is globally asymptotically stable for the reaction-diffusion system 
(II. 2)1 without cross-diffusion. 


3. Effects of cross-diffusion on Turing instability 


For simplicity, we denote 

( 


K(u) = 


(k u + k l3 u 3 )ui 


\ 


(&22 + k2 3 U 3 )u2 

y (k 3 iUi + k 32 U2 + k 33 )u 3 j 


F(u) = 


^ aui(l — Ui) — u±u 3 ^ 
bu 2 ( 1 - u 2 ) - u 2 u 3 
y -cul + (dui + eu 2 )u 3 j 
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Then the reaction-diffusion system (11.21) can be rewritten in matrix notation as: 


<9ii 

— - AK(u) = F(u) 

in 12 x ( 0 , oo), 

9u 

O 

|c 6 

on 12 x ( 0 , oo), 

u(x, 0 ) = (u 10 (x),u 20 (x),u 30 (x)) T 

in 12 . 


(3.1) 


Linearizing the reaction-diffusion system (13.11) about the positive equilibrium u = (u\, u 2 , M 3 ), 
we have 


<9T 

dt 


K u (u)AT = G u (u)T 


(3.2) 


where T = (- 01 , z/> 2 , Vb) T and 


KJn) = 


^ hi + fcl3«3 0 ki3Ul ^ 

0 k 2 2 + ^23^3 ^23^2 

^31^3 ^32^3 &33 + h\Ui + k 32 U 2 ) 


( 


G u (u) 


—alii 0 

0 — bu 2 

du 3 eu 3 


—«2 
-cw 3 y 


Let 0 = fii < /i 2 < /13 < • • • be the eigenvalues of the operator —A on 12 with the 
homogeneous Neumann boundary condition, and E(Hi) be the eigenspace corresponding 
to /it in G 2 (fi). Let X = {u <E [C' 1 (0)] 3 | = 0 on 912}, {0y}j=i l2 ,...,dim£;( w ) be an 

orthonormal basis of E(fa), and = {c0 ij7 -1 c G R 3 }. Then 


oo dim E(ih) 

X = 0X, and Xi= 0 Xij. 

i=l j=1 

For each z > 1, X,; is invariant under the operator K u (u)A + G u (u). Then problem (j3.2[) 
has a non-trivial solution of the form T = c^exp(At) if and only if (A, c) is an eigenpair 
for the matrix — /z,;K u (u) + G u (u), where c is a constant vector. Then the equilibrium u 
is unstable if at least one eigenvalue A has a positive real part for some Hi- 
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The characteristic polynomial of — /XjK u (u) + G u (u) is given by 

Pi (A) = A 3 + A 2i \~ + A u \ + A 0i , (3-3) 

where 

An = (fcn + &22 + fc.33 + ^13^3 + k 23 u 3 + A: 3 i«i + k 32 u 2 )ui + ahi + bu 2 + ch 3 (3.4) 

-^ii — [(&n + ki 3 u 3 )(k 22 + fc 2 3 M 3 + fc 33 + + k 32 u 2 ) + [k 22 + k 23 u 3 ) 

(^33 + k 3 \U\ + k 32 u 2 ) — k 23 k 32 u 2 u 3 — ki 3 k 3 iUiu 3 ]^ 

+ [{bu 2 + cu 3 ) (fcn + ki 3 u 3 ) + au\{k 22 + fc 33 + fc 2 3 h 3 + fc 3 1 u 1 + k 32 u 2 ) 
+ 6 m 2 (^33 + k 3 \U\ + k 32 u 2 ) + cu 3 {k 22 + k 23 u 3 ) + u 2 u 3 {k 23 e — k 32 ) 
+uiu 3 (ki 3 d - k 3 i)\fii + aui(bu 2 + ch 3 ) + bcu 2 u 3 + eh 2 h 3 + duih 3 (3.5) 

^ 0 * — [(^ 33^11 + k 33 k\ 3 u 3 + fc 3 ifcnWi + k 3 \U\U 3 + k 32 knu 2 + k 32 k\ 3 u 2 u 3 )k 22 
A{k 33 kn + k 33 k\ 3 u 3 + fc 3 ifciihi)A: 23 ]' u f + [cm 3 (A;ii + ki 3 u 3 ){k 22 + k 23 u 3 ) 

+ (^33 + ^31 + fc 32 'U 2 )[a'Ui(A :22 + k 23 u 3 ) + bu 2 (ku + fci 3 h 3 )] 

+(&n + k 13 u 3 )(ek 23 u 2 u 3 ) - ak 23 k 32 uiu 2 u 3 + dk 13 ihu 3 (k 22 + fc 23 tx 3 ) 
—ki 3 k 3 ibu\u 2 u 3 — k 3 \k 22 u\u 2 — k 3 ik 32 uiul)ii^ + [ac'Ui'U 3 (fc 2 2 + k 23 u 3 ) 
+bcu 2 u 3 (k u + fci 3 M 3 ) + abuiu 2 {k 33 + k 31 u 3 + k 32 u 2 ) + eu 2 u 3 (kn 
+k l3 u 3 ) + aui(ek 23 u 2 u 3 - k 23 u 2 u 3 ) + kuUiu 2 u 3 bd 

+du 1 u 3 (k 2 2 + k 23 u 3 - fc 3 i6-ui'u 2 h 3 )]-u i + (abc + ae + bd)uiu 2 u 3 . (3.6) 

Let Aij, A 2 j, A 3i be the three roots of (13.31) . In order to obtain the stability of u, we need 
to show that three exists a positive constant S such that 

Re{\u}, Re{\ 2i }, Re{\ 3i } < — <5, for all i> 1. (3.7) 

The aim of the following theorem is to prove that the diffusion alone (with out cross¬ 
diffusion, i.e k 3 1 = k\ 3 = k 32 = k 23 = 0) can not drive instability for this model. 
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Theorem 3.1. Suppose that (12.3p holds and k\ 3 = fc 31 = k 2 3 = k 32 = 0. Then the positive 
equilibrium u 0/ (13.11) is linearly stable. 


Proof: Substituting ki 3 =k 3 i=k 23 =k 32 = 0 into (13.41) . (13.5ft and (13. 6 ft we have 

A- 2 i — + bu 2 + CU 3 + (fcn + fc 22 + &33 )p% > 0 

= (knk22 + knk 33 + k 22 k 33 )p 2 + [a(fc 22 + /c 33 )ui + h(hn + h 33 )h 2 + c(k u + k 2 2)u 3 \pi 
+abu\U2 + acu\u 3 + du 3 u 3 + bcu 2 u 3 + e« 2 «3 > 0 
^ 4 oi = kuk 22 k 33 p 3 + (fcii k 22 cu 3 + knk 33 bu2 + k 22 k 33 aui) p 2 

+(abk 33 u 1 u 2 + dk 22 uiu 2 + ack 22 U\U 3 + bcknu 2 u 3 + eknu 2 u 3 + bdknUiU 2 u 3 ) pi 
+ (ae + bd + abc)uiu 2 u 3 > 0 . 

A direct calculation shows that A 2i A u — A 0 j > 0 for all i > 1. It follows from Routh- 
Hurwize criterion that all the three roots Aij, A 2 %, A 3 j of Pi(X) = 0 have negative real parts 
for each i > 1 . 

Let A = p,i£, then 

Pi(X) = p 3 £ 3 + A 2i p 2 £ 2 + Aup^ + Aq i = Pi(£). 

Since pi —* 00 , as i —* 00 , we have 

p(£) = li m ~ = £ 3 + (fcn + h 22 + k 33 )^ 2 + {k\\k 2 2 + k 22 k 33 + ^11^33)^ + kuk 2 2 k 33 . 

WOO p. 

Applying the Routh-Hurwitz criterion it follows that the three roots £ 1 , £ 2 , £3 of p(£) = 0 all 
have negative real parts. Thus, there exists a positive constant 5 such that Re{£i}, Re{£ 2 }, 
^{£ 3 } < —25. By continuity, we see that there exists io > 1, such that p^ > 1 and the 
three roots £ a , £* 2 , £*3 of p,(£) = 0 satisfy Re{£a}, Re{£, :2 }, Re{£ i3 } < -prf < -p i0 5 < -5 
for any i > i 0l Let —5 = maxi<j<j 0 {i2e{Aii}, Re{Aj 2 }, Re{Aj 3 }} and 5 = min{<5, 5}, Then 
(13.71) holds. Consequently the equilibrium u is linearly stable. 

Note that A 2i > 0, A u > 0, A 0i > 0, and A 2i A u — A 0i > 0 if k 3l = k 32 = 0 since 
the possible negative terms all involve either fc 31 or k 32 . By the same arguments as in 
Theorem 13.11 we have 











Theorem 3.2. Suppose that (j2.3[) holds and k 3l = k 32 = 0, Then the positive equilibrium 
U 0/(0 D is linearly stable. 

Next we consider the Turing instability i.e. the stability of the positive equilibrium 
u = (ui,u 2 ,u 3 ) changing from stable, for the ODE dynamics (11.11) . to unstable for the 
PDE dynamics (11.21) . Here we give sufficient conditions on cross-diffusion which drives 
the instability. fc 31 and k 32 are chosen as variation parameters. 

Theorem 3.3. (1) Suppose that au 3 — u 3 < 0. Consider k 31 as the variation parameter, 
then there exists a positive constant 5 3 i such that when k 33 > S 3 i, the equilibrium u is 
linearly unstable for some domain 0. 

(2) Suppose that bu 2 — u 3 < 0. Consider k 32 as the variation parameter, then there exists 
a positive constant S 32 such that when k 32 > S 32 , the equilibrium u is linearly unstable for 
some domain 0. 

Proof: Denote 


^4(p) ~ — (C 3 /i 3 + C 2 p~' + Ci/a + C 0 ) (3.8) 

where 

C 3 = {(k 33 kn + k 33 ki 3 u 3 + k 3 iknUi + k 3 \UiU 3 + k 32 knu 2 + k 32 ki 3 u 2 u 3 )k 22 
+{k 33 kn + k 33 k\ 3 u 3 + k 3 iknUi)k 23 \ 

C 2 = [cu 3 (k u +ki 3 u 3 )(k 22 +k 23 u 3 ) + (k 33 +k 31 +k 32 u 2 )[au 1 (k 22 +k 23 u 3 ) 
+bu 2 (k u + h 3 u 3 )\ + (k u + k 13 u 3 )(ek 23 u 2 u 3 ) - ak 23 k 32 uiu 2 u 3 
+dk l3 u 1 u 3 {k 22 + k 23 u 3 ) - fci 3 /c 3 i6hih 2 h 3 - k 31 k 22 uiu 2 - k^k^ul] 

Ci = [abuiu 2 k 31 (aui - u 3 ) + auiu 2 k 32 {bu 2 - u 3 ) + acuiu 3 (k 22 + k 23 u 3 ) + 
bcu 2 u 3 (k u + k 13 u 3 ) + abuiu 2 k 33 + eu 2 u 3 (k u + k 13 u 3 ) + aek 23 uiu 2 u 3 
+knbduiu 2 u 3 + duiu 3 (k 22 + k 23 u 3 ) 

C 0 = (abc + ae + bd)uiu 2 u 3 . 
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Case 1: k 3 1 is the variation parameter. 

We assume that au\ — u 3 < 0. The following arguments by continuation are based on 
the fact that each root of the algebraic equation (13.8ft is a continuous function of the 
variation parameter k 3 \. It is easy to prove that equation (13.81) has three real roots 
Hi ' = /i| ?) (/c 3 i), i — 1,2,3 when fc 31 goes to infinity and they are lim^^oo p^\k 3 i) < 0, 
linifcg^oo /i| 2) (h 3 i) = 0 and lim^^oo pf\k 3 i) > 0. By continuation, there exists a positive 
constant <5 3 i such that when k 33 > d 3 i, C\ > 0 and det(H(/i)) has three real roots. Because 
C 3 > 0 and C 0 > 0, the mumber of sign changes of (j3.8[) is exactly two. Therefore by 
Descartes’rule, the three real roots have the following properties: 

( 1 ) — oo < p^ < 0 < p^ < pf^ < oo, 

(2) det(y4(/i)) > 0 if p G (-oo,^) U (PiKp?), 

(3) det(y4(p)) < 0 if p G {p^\ pf^) U (pf\ oo). 

If /ij G (pf\pf^) for some z, then det(H(/Zj)) > 0 by (2), and consequently A 0i = 
— det(Ap)) < 0. The number of sign of changes of the characteristic polynomial (13.31) 
Pi{ A) = X 3, + A 2 i\ 2 + Au\ + Aoi is either one or three. By Descartes’rule, the characteristic 
polynomial (13.31) has at least one positive eigenvalue. Hence, the equilibrium u of (11.21) 
is linearly unstable for any domain fl on which at least one eigenvalue pi of —A is in the 
interval {pf\pf^). 

Case 2 : k 32 is the variation parameter. 

We assume that bu 2 — zz 3 < 0. The following arguments by continuation are based on 
the fact that each root of the equation (j3.8[) is a continuous function of the variation k 32 . 
It is easy to prove that equation (j3.8[) has three real roots p^ = (^. 32 ), i = 1,2,3 

when k 32 goes to infinity and they are lim^^oo p^\k 32 ) < 0 , lim^^,*, p^\k 32 ) = 0 and 
linifcg^oo/4 7 (k 32 ) > 0. By continuation, there exists a positive constant <5 32 such that 
when k 32 > <5 32 , C 3 > 0 and det(H(/i)) has three real roots. Because C 3 > 0 and Co > 0, 
the mumber of sign changes of (13.8j) is exactly two. Therefore by Descartes’rule,the three 
real roots have the folloeing properties: 

( 1 ) — 00 < p ^ < 0 < p ^ < p ^ < 00, 

(2) det {A(p)) > 0 if p G (- 00 , p% ] ) \J(p 2 \ p 2 ] ), 
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(3) det (A(fi)) < 0 if /i G (ix 2 ] , U(/4 3) »°°)- 

If i-ii G (/4 2 \a 4^) f° r some i , then det(v4(/Zj)) > 0, and consequently A 0i = — det(A(q) < 
0. By similar argument as case 1, The number of sign of changes of the characteristic 
polynomial (13.3|) pfX) = A 3 + A 2 iX 2 + A U X + A 0i is either one or three. By Descartes’rule, 
the characteristic polynomial (j3.3[) has at least one positive eigenvalue. Hence, the equi¬ 
librium u of (11.2ft is linearly unstable for any domain H on which at least one eigenvalue 

/o\ / o\ 

Pi of —A is in the interval (jf 2 ,p 2 ). 

Remark 3.1. (i) In Theorem 3.3, the condition au\ — ^3 < 0 and bu ,2 — u 3 < 0 are 
compatible with the condition (2.3) respectively. 

(ii) k 3 1 and k 32 can be chosen as variation parameters because the number of sign 
of change for the polynomial (3.8) could be bigger than one for large values of fc 31 or 
k 3 2 > By descartes’ rule, the polynomial (3.8) could have positive roots which lead to linear 
instability. 

(in) Biological interpretation: In our model, the third species preys on the first and 
second one. The positive steady state of the model can be broken by the reaction-diffusion 
among two species on the model. Case one: In this case, the first species are assumed to 
reproduce exponentially unless subject to intra-species competitions and predation. This 
exponential growth is represented in the equation by the term au 1 . The level of intra¬ 
species competitions among the first species is assumed to be proportional to the population 
density of first species by the term au\. The rate of predation upon the prey is assumed 
to be proportional to the rate at which the predators and the prey meet by the term u 3 u 3 , 
when the effect on first species due to the fact the third species preys on the first one 
u 3 are larger than the effects on first species due to the intra-species competitions among 
first species au\, the large cross-diffusion of the third species due to the first species k- M 
can break the stability of the positive steady state. In other words, if the predator has a 
dominate effect on the decreasing of the prey such as predation rate is lager than the rate 
of intra-species competitions, then the predator with large cross-diffusion can destabilize 
the constant steady state. Case two: hi this case, the third species shall have a dominate 
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effect on the decreasing of the second species. Because bu 2 — h 3 < 0 implies bu 2 < u 3 , 
predation rate of third species on the second species is large than the rate of intra-species 
competitions in second species. The similar situation as in case one happens in the case 
two: the predator with large cross-diffusion can destabilize the constant steady state. 

4. Numerical calculations 

In this section, using numerical methods, we illustrate that the cross-diffusion induces 
spatial patterns. The initial data is taken as a uniformly distributed random perturbation 
around the equilibrium state (ui,u 2 ,uf) in hi, with a variance lower than the amplitude 
of the final patterns. More precisely, 

u w (x) = ui + 771 (x), u 20 (x ) = u 2 + rj 2 (x), u 30 (x) = u 3 + rj 3 (x), 

where 7 p G [—1.5,1.5] for i = 1,2,3. I 11 view of Theorems 13.11 and 13.31 the Turing 
parameter space is (j2.3[) under which spatial patterns can occur. Thus, in the system 
(11.211 we fix a = 1, b = 1, c = 0.1, d = 0.1, e = 0.1, ku = 0.1, k± 3 = 0.1, k 22 = 0.1, 
k 23 = 0.1, k 3 \ = 0.1, k 33 = 0.1. 

In Figure [U we show that the real part of the eigenvalues A as a function of the cross¬ 
diffusion coefficient d 32 . From the characteristic polynomial of (I3.3f) . we can determine the 
value of d 32 such that Re( A) > 0. Now we will implement some numerical simulations for 
the system (11.211 . The domain is confined to a square domain hi = [0, L x ] x [0 ,L y \ C R 2 . 
The wavenumber for this two dimensional domain is thereby 

k = 7t(t n/L x , n/L y ), and |k| = fm/L x ) 2 + ( n/L y ) 2 , m, n — 0, 1,.... 

We consider system (11.211 in a fixed domain L x = 40 and L y = 40, and resolve it on a grid 
with 100 x 100 sites with the space step of Ax = Ay = 1. For the evolution in time, we 
apply a first order backward Euler time advancing scheme with a time step At = 0.005. 
By discretizing the Laplacian in the grid with lattice sites denoted by (i, j), the nine-point 
formula is 

Att| (ij) = g^[4 ai(i,j)u(i - 1 ,j) +4 a r (i,j)u(i + l,j) +4 a d (i,j)u(i,j - 1) 
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0.05 



Figure 1: Dispersion relation for the real part of the eigenvalues, Re (A) versus the cross-diffusion 
coefficient /C 32 . 
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Figure 2: Bifurcation diagram for Turing onset. Maximum and minimum u\ for different cross-diffusion 
in the transition from the homogeneous state to Turing pattern. 
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Figure 3: Spatial patterns change quantitatively with different k 32 as 1.7, 1.8, 1.9, and 2. The other 
parameters are stated in the text. The steps of the iteration for time is 40000. 

+4 a u (i,j)u(i,j + 1) + ai(i,j)u(i - 1, j + 1) + a u (i, j)u(i + l,j + 1) 
+a d (i,j)u(i - 1, j - 1) + a r (i, j)u(i + 1, j - 1) - 20 u(i,j)\, 

where the matrix elements of ai, a r , ad, a u are unity except at the boundary. When (i,j) is 
at the left boundary, that is i — 0, we define ai(i, j)u(i— 1, j) = u(i+l,j), which guarantees 
zero-flux of reactants in the left boundary. Similarly we define a r (i,j), a,d(i,j), a u (i,j ) 
such that the boundary is no-flux. The nine-point formula for the Laplacian can have a 
one-step error of 0(Ax 4 ). 

In Figure [21 we compare the density of u\ before and after the onset of Turing patterns. 
Results are qualitatively similar for w 2 and u 3 , and hence omitted. In the case of k 32 less 
than 1.6, i.e. the Turing instability does not occur, we see that the density of U\ is 
homogeneous. In the case of k 32 larger than 1.6, i.e. Turing instability happen, we see 
that the density of U\ is spatial inhomogeneous. 

Now we study the change of the spatial patterns qualitatively and quantitatively with 
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different k 32 . In general, the selection of stripe pattern or spot pattern depends upon the 
non-linearities of the reaction kinetics. Specifically, it has been shown that the presence 


of quadratic nonlinearities in the reaction 
of quadratic terms leads to stripe pattern 


dnetics leads to spot pattern, but the absence 
4], Noticing that the reaction kinetics of (II.2p 


only has quadratic nonlinearities, in view of the theory of pattern selection Q, all the 
spatial patterns are spot patterns. In Figure El we also illustrate the quantitative change 
of the spatial patterns with the different k 32 . From this simulations, we can conclude that 
with the increasing of k 32 , the spatial patterns converge to regular spotted patterns. The 
striped patterns can not occur in our model. 


5. Comparisons and conclusions 

In this paper, we have developed a theoretical framework for studying the phenomenon 
of pattern formation in a two-prey one-predator system. Applying a stability analysis and 
suitable numerical simulations, we investigate the Turing parameter space, the associated 
pattern type and the Turing bifurcation diagram. The proposed approach has applicability 
to other reaction-diffusion systems including cross-diffusion, such as chemotaxis and cell 
motility models. In this context, it is of great interest to us the development of a general 
mathematical and numerical framework that allows for the treatment of certain degenerate 
quasilinear parabolic systems modeling bacterial growth, that are known to involve several 
important phenomena such as fractal morphogenesis and branching patterns. 

It is worth mentioning that the authors have also the rokyof cross-diffusion on pattern 


formation for Lotka-Volterra type models in 5), 21]. In 21], by considering a Holling 


Tanner predator-prey model the authors investigated the Turing bifurcation and obtained 


the pattern selection mechanism. In jsj], by studying the h 


attained the spiral patterns. Apart from these work 


opf bifurcation the authors 


2 lj, what our model consider 


is a three species model. The difficulty is that the characteristic equation of our model 
is a cubic equation. We use the continuity of the cubic functions to overcome it. Our 
novelty is that we have obtained the bifurcation diagram for Turing onset by numerical 
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simulations, which shows the transition from the homogeneous steady state to the Turing 
patterns. 

It is well-known that for a classical competitive model, the formation of patterns does 
not occur. We introduce the cross-diffusion into the particular two-prey one-predator 
model, and show that this gives rise to Turing-like spatial patterns. All this is confirmed 
with the help of illustrating numerical simulations. 
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